home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / VideoToolbox / Utilities / Quick3 / Quick3.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-10-13  |  14.9 KB  |  370 lines  |  [TEXT/KAHL]
  1. /*
  2. Quick3.c 
  3. Copyright (c) 1990-1993 Denis G. Pelli 
  4. Quick3.c is a C program that makes maximum likelihood fits of simple models to
  5. psychometric data. This is a replacement for, and to some extent is derived from
  6. the old QUICK.FOR program written by A.B.Watson. See A.B.Watson (1979)
  7. Probability summation over time. Vision Research 19, 515-522.
  8.  
  9. The principal difference between this program and the original QUICK, is in how
  10. the goodness of fit is calculated. Both programs assess the goodness of fit by
  11. comparing the log likelihoods of the model fit (usually Weibull) with a null
  12. hypothesis. Minus two times this log likelihood ratio is approximately chi
  13. square, with a number of degrees of freedom equal to the difference in degrees
  14. of freedom of the two hypotheses. QUICK used an unconstrained psychometric
  15. function, which adopted the actual proportion correct at each contrast. Quick3
  16. uses a monotone hypothesis which is constrained to be monotonically increasing
  17. as a function of contrast. The virtue of the latter hypothesis is that it is
  18. easier to estimate its degrees of freedom when there are unequal numbers of data
  19. points at each contrast. (However, this estimate is seat-of-my-pants.) Finally,
  20. this program actually computes the significance for you, from the chi square
  21. value and the estimate of its degrees of freedom. This significance value is of
  22. some use in choosing a small fraction of your fits to discard (e.g. if
  23. significance is less than 0.05). This significance is misleadingly low when both
  24. hypotheses have nearly the same number of degrees of freedom (i.e. you only have
  25. a few contrasts). Looking at the significance is helpful, but not a substitute
  26. for thinking about whether your results are really reasonable.
  27.  
  28. Note that Quick3.c is just a front end, dealing with the user interface and
  29. reading and writing files. All the work is done by one call to
  30. PsychometricFit(). In many cases it will be convenient to call that routine
  31. directly from your program that collects (or generates) the data, rather than
  32. saving the data in a file for subsequent analysis by Quick3.c. If you want to
  33. use a psychometric model other than the Weibull function all you need to do is
  34. write a small function to implement it, based on Weibull.c. It should be easy.
  35. PsychometricFit() uses whatever psychometric function is supplied in the call.
  36.  
  37. Quick3 reads in your data from a text file and analyzes them. The results are
  38. shown on the screen, and are saved in two kinds of output file. The .fit file is
  39. in Excel format and has just a minimal one-line summary for each condition,
  40. giving the parameters and goodness of fit. The .plot file is in text format and
  41. is suitable for plotting (e.g. by CricketGraph) of the psychometric data along
  42. with the Weibull and Monotonic fits. For Macintosh users, the output files are
  43. TEXT files, but have the appropriate creator so that double-clicking them opens
  44. the data file in the appropriate application, either Excel or Cricket Graph.
  45. For maximum convenience, change the definition of PLOT_CREATOR to match that of 
  46. your favorite graphing program.
  47.  
  48. You should also have received a CricketGraph format files called "CricketGraph
  49. Quick3" (for the original CricketGraph) and "CA-CricketGraphIII Quick3" (for the new
  50. CricketGraphIII). Put the appropriate file in the same folder as the
  51. CricketGraph application. When CricketGraph starts up it will read in the format
  52. file, and later you'll be able to select that format from CricketGraph's Format
  53. menu when you plot your data.
  54.  
  55. Quick3's input file format is simple and flexible. It is a text file. Lines
  56. beginning with "#" are treated as comments. Comments are either copied directly
  57. to the .fit file, or used as a name for the next "condition" (i.e. block of
  58. data). The file must begin with at least one comment line. The last of several
  59. contiguous comment lines will be taken to be the name for the following
  60. condition. All the non-comment lines until the next comment will be interpreted
  61. as the data for that condition.
  62.  
  63. An experiment usually consists of several conditions, or blocks of trials.
  64. Within the condition the observer will be tested repeatedly at various
  65. contrasts. The results consist of the number of trials correct at each contrast.
  66. The results may be described by a set of "contrast records". Each contrast
  67. record consists of a contrast (e.g. 0.012), the number of trials at that
  68. contrast (e.g. 1 or 100), and the number of correct responses (e.g. 0 or 79).
  69. Quick3 expects a contrast record to look like this:
  70. 0.012 100 79
  71. The three numbers are all on one line, and separated by white space, so sscanf
  72. can parse them. Anything on the line after the third number is ignored. All
  73. contrast records until the next comment (line beginning with "#") or the end of
  74. file are assumed to be from the same condition and are analyzed together. The
  75. contrast records may be in any order. The records will be sorted into order of
  76. ascending contrast. You may have any number of trials at each contrast, up to
  77. about two billion (i.e. LONG_MAX in limits.h). It's okay to have multiple
  78. records with equal contrasts. The data will be merged, adding up the trials at
  79. that contrast. You can even be really lazy and just write out each trial as a
  80. separate contrast record, in whatever order, without keeping track of how many,
  81. etc.
  82.  
  83. You may add records with zero trials (and zero correct) so as to have fitted
  84. values computed at those contrasts in the .plot file.
  85.  
  86. Here's an example of a data file, suitable for analysis by Quick3:
  87. #This is a sample file. 
  88. #The first condition is called "monocular" and has 60 trials at 6 contrasts.
  89. #monocular
  90. 0.1 10 4
  91. 0.2 10 5
  92. 0.3 10 7
  93. 0.4 10 6 Anything after the specified data is ignored.
  94. 0.5 10 8
  95. 0.6 10 10
  96. # The "binocular" condition is bigger: 1000 trials.
  97. #binocular
  98.  0.056     100    53
  99.  0.064     100    53
  100.  0.073     100    65
  101.  0.083     100    77
  102.  0.094     100    81
  103.  0.107     100    84
  104.  0.121     100    96
  105.  0.138     100    99
  106.  0.156     100   100
  107.  0.178     100   100
  108. #all done
  109.  
  110. Here's the resulting .fit file (tabs will space properly in Excel, but not here):
  111. Condition    logAlpha    beta    gamma    delta    free params    signif.    Chi sq.    d.f.    trials    contrasts
  112. #This is a sample file. 
  113. #The first condition is called "monocular" and contains 60 trials at 6 contrasts.
  114. monocular    -0.305    8.269    0.524    0    4    0.183    1.773    1    60    6
  115. # The "binocular" condition is bigger: 1000 trials.
  116. binocular    -1.015    3.811    0.482    0    4    0.479    5.521    6    1000    10
  117. #all done                                        
  118.  
  119. SOURCES:
  120.  
  121. Quick3.h
  122. LogLikelihood.c
  123. MonotonicFit.c
  124. PsychometricFit.c
  125. Quick3.c
  126. SortAndMergeContrasts.c
  127. Weibull.c
  128. #From Denis Pelli's VideoToolbox:
  129. VideoToolbox.h
  130. Binomial.c
  131. ChiSquare.c
  132. Normal.c
  133. SetFileInfo.c        # Used only on the Macintosh
  134. #From Numerical Recipes in C:
  135. nr.h
  136. NRUTIL.h
  137. BRENT.C
  138. F1DIM.C
  139. LINMIN.C
  140. MNBRAK.C
  141. NRUTIL.C
  142. POWELL.C
  143.  
  144. HISTORY:
  145. 4/7/90        dgp    wrote it.
  146. 4/10/90        dgp    changed .fit format from %.3f to %.4f. Accept only printing characters
  147.                 in the condition name, ignoring everything after the first non-printing
  148.                 character. This deals with the fact that if an Excel file is used for
  149.                 input, there are lots of trailing tabs that should be ignored. Checked
  150.                 that files were actually open before closing 'em. Jeesh!
  151. 10/29/90    dgp    tidied up comments.
  152. 1/20/91        dgp    updated calls to BinomialUpperBound & BinomialUpperBound to
  153.                 conform to new definition.
  154. 8/24/91        dgp    Made compatible with THINK C 5.0.
  155. 11/17/92    dgp Now SetVol() after calling SFGetFile, so that file can be in any folder,
  156.                 not just the same folder that Quick3 is in.
  157. 1/19/93        dgp    put #ifs around the Mac-dependent code.
  158. 2/20/93        dgp    added call to Require().
  159. 8/5/93        dgp    moved call to Require() into a separate main(), since application
  160.                 failed on Mac Plus before it got to the call to Require(). Thanks
  161.                 to Robert Friedman friedman@cortex.health.ufl.edu for reporting the
  162.                 problem.
  163. 10/13/93    dgp    added comments about CricketGraphIII. Allow compile-time choice of
  164.                 Mac file's creator.
  165. */
  166. //#define PLOT_CREATOR    'CGRF'    /* Create Macintosh Cricket Graph file */
  167. #define PLOT_CREATOR    'CRGR'    /* Create Macintosh CA-Cricket Graph III file */
  168. #include "VideoToolbox.h"
  169. #include "Quick3.h"
  170. #include <nr.h>                /* Numerical Recipes in C*/
  171. #if MACINTOSH
  172.     #include <QuickDraw.h>
  173.     #include <StandardFile.h>
  174. #endif
  175. char *MyFgets(char *s, int length, FILE *stream);
  176. void Quick3(void);
  177.  
  178. void main(void)
  179. {
  180.     #if MACINTOSH
  181.         Require(0);    // check for any required cpu and fpu.
  182.     #endif
  183.     Quick3();
  184. }
  185. void Quick3(void)
  186. {
  187.     static dataRecord data,monotonicData;
  188.     contrastRecord *cPtr;
  189.     static paramRecord params, initParams;
  190.     double *paramPtr=¶ms.logAlpha;    /* a generic way of accessing the parameters */
  191.     int i;
  192.     double chiSquare,modelLL,monotonicLL;    /* log likelihood */
  193.     int chiSquareDF,modelDF,monotonicDF;    /* degrees of freedom */
  194.     double significance;
  195.     FILE *dataFile=NULL,*fitFile=NULL,*plotFile=NULL;
  196.     static char dataFileName[50],fitFileName[50],plotFileName[50],string[100];
  197.     static char conditionName[50];
  198.     long trials;
  199.     unsigned int a;
  200.     PsychometricFunctionPtr ModelFunction=&Weibull;    /* a psychometric function */
  201.     static const char modelName[]="Weibull";
  202.     static const char paramName[PARAMS][16]={"logAlpha","beta","gamma","delta"};
  203.     #if MACINTOSH
  204.         Point where;
  205.         static SFTypeList typeList;
  206.         static SFReply reply;
  207.     #endif
  208.  
  209.     /* initial values for the parameters of the psychometric function */
  210.     params.logAlpha=0.0;
  211.     params.beta=3.5;
  212.     params.gamma=0.5;
  213.     params.delta=0.01;
  214.     printf("\n");    /* ask THINK C to init QuickDraw */
  215.     
  216.     #if MACINTOSH
  217.         where.h=100;
  218.         where.v=100;
  219.         typeList[0]='TEXT';
  220.         reply.version=0;
  221.         SFGetFile(where,"\p",NULL,1,typeList,NULL,&reply);
  222.         if(!reply.good)PrintfExit("Couldn't open file.\n");
  223.         SetVol(NULL,reply.vRefNum);    /* look in that folder */
  224.         sprintf(dataFileName,"%#s",reply.fName);
  225.     #else
  226.         // Insert code here to get data file name into dataFileName.
  227.     #endif
  228.     dataFile=fopen(dataFileName,"r");
  229.     if(dataFile==NULL) {
  230.         PrintfExit("\007Sorry, can't find file \"%s\".\007\n",dataFileName);
  231.     }
  232.     printf("Reading \"%s\"\n",dataFileName);
  233.     strcpy(string,dataFileName);
  234.     if(strstr(string,".data"))strcpy(strstr(string,".data"),"");
  235.     else if(strstr(string,".yes"))strcpy(strstr(string,".yes"),"");
  236.     if(strlen(string)>0){
  237.         sprintf(fitFileName,"%s.fit",string);
  238.         printf("Creating output files:\n%s\n%s.<condition name>.plot\n",fitFileName,string);
  239.         fitFile=fopen(fitFileName,"w");
  240.         if(fitFile==NULL)PrintfExit("Sorry, I can't create file \"%s\".\007\n",fitFileName);
  241.         #if MACINTOSH
  242.             SetFileInfo(fitFileName,'TEXT','XCEL');        /* Excel file */
  243.         #endif
  244.     }
  245.     modelDF=PARAMS;    /* number of parameters to be adjusted in fitting */
  246.     printf("How many parameters shall be adjustable? 0 to %d (%d):",PARAMS,modelDF);
  247.     gets(string);
  248.     sscanf(string,"%d",&modelDF);
  249.     for(i=1;i<modelDF;i++){
  250.         printf("Initial value for %s? (%.2f):",paramName[i],paramPtr[i]);
  251.         gets(string);
  252.         sscanf(string,"%lf",¶mPtr[i]);
  253.     }
  254.     for(i=modelDF;i<PARAMS;i++){
  255.         printf("Fixed value for %s? (%.2f):",paramName[i],paramPtr[i]);
  256.         gets(string);
  257.         sscanf(string,"%lf",¶mPtr[i]);
  258.     }
  259.     initParams=params;
  260.     
  261.     printf("Reading %s\n",dataFileName);
  262.     if(fitFile)fprintf(fitFile,"Condition\t"
  263.         "%s\t%s\t%s\t%s\tfree params"
  264.         "\tsignif.\tChi sq.\td.f.\ttrials\tcontrasts\n",
  265.         paramName[0],paramName[1],paramName[2],paramName[3]);
  266.     while(!feof(dataFile)){
  267.         MyFgets(string,sizeof(string),dataFile);
  268.         printf("\n");
  269.         while (a=fgetc(dataFile),ungetc(a,dataFile),a=='#' || a==EOF){
  270.             printf("%s\n",string);                                /* echo comment */
  271.             if(fitFile)fprintf(fitFile,"%s\n",string);
  272.             MyFgets(string,sizeof(string),dataFile);
  273.             if(feof(dataFile))break;
  274.         }
  275.         /* get condition name. Strip leading # and trailing nonprinting junk */
  276.         strcpy(conditionName,&string[1]);
  277.         for(i=0;i<strlen(conditionName);i++){
  278.             if(!isprint(conditionName[i]))conditionName[i]=0;
  279.         }
  280.         printf("\n");
  281.         data.contrasts=i=0;
  282.         while(a=fgetc(dataFile),ungetc(a,dataFile),a!='#' && a!=EOF){    /* read data */
  283.             if(i==MAX_CONTRASTS){
  284.                 SortAndMergeContrasts(&data);
  285.                 if(i<MAX_CONTRASTS)break;
  286.                 PrintfExit("Quick3: Sorry, too many different contrasts. >MAX_CONTRASTS==%d\007\n",
  287.                     MAX_CONTRASTS);
  288.             }
  289.             MyFgets(string,sizeof(string),dataFile);
  290.             sscanf(string,"%lf\t%ld\t%ld",
  291.                 &data.c[i].contrast,&data.c[i].trials,&data.c[i].correct);
  292.             i++;
  293.             data.contrasts=i;
  294.         }
  295.         if(data.contrasts==0)break;
  296.         SortAndMergeContrasts(&data);
  297.         trials=0;
  298.         for(i=0;i<data.contrasts;i++)trials+=data.c[i].trials;    /* total trials */
  299.         params=initParams;    /* initial guess & fixed values */
  300.         if(modelDF>0)
  301.             params.logAlpha=log(data.c[data.contrasts/2].contrast)/log(10.0);    /* median contrast */
  302.         printf("%s\n",conditionName);
  303.         printf("\t%s   %s   %s   %s  contrasts trials signif. Chi sq.   d.f.\n",
  304.             paramName[0],paramName[1],paramName[2],paramName[3]);
  305.         printf("Guess:");
  306.         printf("\t%7.2f%8.1f%8.2f%8.2f\n",
  307.             paramPtr[0],paramPtr[1],paramPtr[2],paramPtr[3]);
  308.         significance=PsychometricFit(¶ms,ModelFunction,&data,&modelLL,modelDF,
  309.             &chiSquare,&chiSquareDF);
  310.         printf("Fit:");
  311.         printf("\t%7.2f%8.1f%8.2f%8.2f",
  312.             paramPtr[0],paramPtr[1],paramPtr[2],paramPtr[3]);
  313.         printf("%8d%8ld",data.contrasts,trials);
  314.         printf("%8.2f%8.1f%8d\n",significance,chiSquare,chiSquareDF);
  315.         if(fitFile)fprintf(fitFile,"%s"
  316.                 "\t%7.4f\t%7.4f\t%7.4f\t%7.4f\t%7d"
  317.                 "\t%7.4f\t%7.4f\t%7d"
  318.                 "\t%7ld\t%7d\n",
  319.                 conditionName,
  320.                 paramPtr[0],paramPtr[1],paramPtr[2],paramPtr[3],modelDF,
  321.                 significance,chiSquare,chiSquareDF,
  322.                 trials,data.contrasts);
  323.     
  324.         /* Now create plot file */
  325.         if(fitFile){
  326.             strcpy(plotFileName,fitFileName);
  327.             plotFileName[strlen(plotFileName)-strlen(".fit")]=0;    /* strip off the ".fit" */
  328.             sprintf(plotFileName,"%s.%s.plot",plotFileName,conditionName);
  329.             plotFile=fopen(plotFileName,"w");
  330.             if(plotFile==NULL){
  331.                 PrintfExit("Sorry, I can't create file \"%s\".\n\007",plotFileName);
  332.             }
  333.             #if MACINTOSH
  334.                 SetFileInfo(plotFileName,'TEXT',PLOT_CREATOR);    /* specify graphing program */
  335.                 if(PLOT_CREATOR=='CGRF')fprintf(plotFile,"*\n");
  336.             #endif
  337.             monotonicData=data;
  338.             MonotonicFit(&monotonicData,&monotonicLL,&monotonicDF);    /* overwrites data with fit */
  339.             fprintf(plotFile,"Contrast\t%s\tLower bound\tUpper bound\t%s fit\tMonotone fit\tTrials\tCorrect\tParameters\n",
  340.                 conditionName,modelName);
  341.             for(i=0;i<data.contrasts;i++){
  342.                 cPtr=&data.c[i];
  343.                 fprintf(plotFile,"%.3f",cPtr->contrast);
  344.                 if(cPtr->trials>0){
  345.                     fprintf(plotFile,"\t%.3f",cPtr->correct/(double)cPtr->trials);
  346.                     fprintf(plotFile,"\t%.3f\t%.3f",
  347.                         BinomialLowerBound(0.95,cPtr->correct,cPtr->trials),
  348.                         BinomialUpperBound(0.95,cPtr->correct,cPtr->trials));
  349.                 }
  350.                 else
  351.                     fprintf(plotFile,"\t\t\t");
  352.                 fprintf(plotFile,"\t%.3f",(*ModelFunction)(cPtr->contrast,¶ms));
  353.                 if(cPtr->trials>0)
  354.                     fprintf(plotFile,"\t%.3f\t%5ld\t%5ld",
  355.                         monotonicData.c[i].correct/(double)monotonicData.c[i].trials,
  356.                         cPtr->trials,cPtr->correct);
  357.                 else fprintf(plotFile,"\t\t\t");
  358.                 if(i<PARAMS)fprintf(plotFile,"\t%s=%.3f",paramName[i],paramPtr[i]);
  359.                 if(i==0)fprintf(plotFile,"\talpha=%.4f",pow(10.0,paramPtr[0]));
  360.                 if(i==PARAMS)fprintf(plotFile,"\tsignificance=%.3f",significance);
  361.                 fprintf(plotFile,"\n");
  362.             }
  363.             if(plotFile)fclose(plotFile);
  364.         }
  365.     }
  366.     if(fitFile)fclose(fitFile);
  367.     if(dataFile)fclose(dataFile);
  368. }
  369.  
  370.